Step 1. Setting Up the Characteristic Equations

For a first–order PDE of the form \(\displaystyle P(x,y,z)z_x + Q(x,y,z)z_y = R(x,y,z) \), the characteristic equations are \begin{align} \frac{dx}{ds} &= P(x,y,z), \tag{1} \\ \frac{dy}{ds} &= Q(x,y,z), \tag{2} \\ \frac{dz}{ds} &= R(x,y,z). \tag{3} \end{align} For our PDE, we identify \begin{align} P(x,y,z) &= x(z+2a), \tag{4} \\ Q(x,y,z) &= xz+2yz+2ay, \tag{5} \\ R(x,y,z) &= z(z+a). \tag{6} \end{align} Hence, the characteristic system is \begin{align} \frac{dx}{ds} &= x(z+2a), \tag{7} \\ \frac{dy}{ds} &= xz+2yz+2ay, \tag{8} \\ \frac{dz}{ds} &= z(z+a). \tag{9} \end{align}

Step 2. First Invariant: Eliminating s Between x and z

We eliminate the parameter \( s \) by writing \[ \frac{dx}{dz} = \frac{\frac{dx}{ds}}{\frac{dz}{ds}} = \frac{x(z+2a)}{z(z+a)}. \] Separating variables gives \[ \frac{dx}{x} = \frac{z+2a}{z(z+a)}\,dz. \]

To integrate the right-hand side, we use partial fractions. Assume \[ \frac{z+2a}{z(z+a)} = \frac{A}{z} + \frac{B}{z+a}. \] Multiplying by \(z(z+a)\) yields \begin{align} z+2a &= A(z+a) + Bz \nonumber\\[6pt] &= (A+B)z + Aa. \tag{10} \end{align} Equating coefficients we have \begin{align} A+B &= 1, \tag{11} \\ Aa &= 2a \quad\Longrightarrow\quad A = 2 \quad (\text{assuming } a\neq0). \tag{12} \end{align} Thus, \( B = 1-2 = -1 \) and \[ \frac{z+2a}{z(z+a)} = \frac{2}{z} – \frac{1}{z+a}. \]

Integrate both sides: \begin{align} \int \frac{dx}{x} &= \int \left(\frac{2}{z} – \frac{1}{z+a}\right)dz \nonumber\\[6pt] \ln|x| &= 2\ln|z| – \ln|z+a| + C_1. \tag{13} \end{align} Exponentiating, we get \[ x = K_1\,\frac{z^2}{z+a}, \quad \text{with } K_1 = e^{C_1}. \] Rearranging leads to the first invariant: \[ \frac{x(z+a)}{z^2} = \text{constant}. \tag{14} \]

Step 3. Second Invariant: Eliminating s Between y and z

Next, consider \[ \frac{dy}{dz} = \frac{\frac{dy}{ds}}{\frac{dz}{ds}} = \frac{xz+2yz+2ay}{z(z+a)}. \] Notice that \[ xz+2yz+2ay = xz + 2y(z+a), \] so we can write \[ \frac{dy}{dz} = \frac{x}{z+a} + \frac{2y}{z}. \tag{15} \]

Substitute the expression for \( x \) from (13): \[ x = K_1\,\frac{z^2}{z+a} \quad\Longrightarrow\quad \frac{x}{z+a} = K_1\,\frac{z^2}{(z+a)^2}. \] Thus, equation (15) becomes \[ \frac{dy}{dz} = K_1\,\frac{z^2}{(z+a)^2} + \frac{2y}{z}. \tag{16} \] Write this linear ODE in standard form: \[ \frac{dy}{dz} – \frac{2}{z}\,y = K_1\,\frac{z^2}{(z+a)^2}. \tag{17} \]

The integrating factor is \[ \mu(z) = \exp\Bigl(-\int \frac{2}{z}\,dz\Bigr)=z^{-2}. \tag{18} \] Multiplying (17) by \( z^{-2} \) gives \[ \frac{d}{dz}\Bigl(y\,z^{-2}\Bigr) = \frac{K_1}{(z+a)^2}. \tag{19} \] Integrate with respect to \( z \): \begin{align} y\,z^{-2} &= K_1\int \frac{dz}{(z+a)^2} + C_2 \nonumber\\[6pt] &= -\frac{K_1}{z+a} + C_2, \tag{20} \end{align} where \( C_2 \) is an arbitrary constant.

Multiplying (20) by \( z^2 \) results in \[ y = -\frac{K_1\,z^2}{z+a} + C_2\,z^2. \tag{21} \] But since \[ x = K_1\,\frac{z^2}{z+a}, \] we have \[ -\frac{K_1\,z^2}{z+a} = -x. \] Therefore, \[ y = -x + C_2\,z^2 \quad\Longrightarrow\quad \frac{x+y}{z^2} = C_2. \tag{22} \] This is the second invariant.

Step 4. The Integral Surface

Since the invariants remain constant along characteristics, the general integral surface is given by an arbitrary function of these invariants: \[ \boxed{\Phi\!\Biggl(\frac{x(z+a)}{z^2},\, \frac{x+y}{z^2}\Biggr) = 0.} \tag{23} \]